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ABSTRACT 

Signal analysts have traditionally relied on the Discrete 
Fourier Transform and various data windowing schemes for 
Signal detection and classification. Some signals, notably 
those of a transient nature, are inherently difficult to 
analyze with these traditional tools. The Discrete Wavelet 
Transform has recently generated considerable interest in 
several areas of digital signal processing and a determination 
of its suitability as a signal analysis tool is necessary. 
Associated with wavelet theory is the concept’ of 
multiresolutional analyses which allow examination of a signal 
at different scales. 

This thesis investigates dyadic discrete wavelet 
decompositions of signals. A new multiphase wavelet transform 
is proposed and investigated. The - multiphase transform 
technique is shown to be useful in transient signal analysis. 
Several MATLAB” programs that perform multiresolutional 


analyses with various supporting features are provided. 
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I. INTRODUCTION 


The Discrete Fourier Transform has traditionally been the 
dominant tool in digital signal processing for extracting 
waveform features required in the analyses of signals. 
Characteristic of the Fourier transform is the global manner 
in which it operates on the input signal, thereby rendering 
this technique incapable of providing time localization of 
frequency components. The analyses of signals that are 
transitory in nature will particularly suffer the deleterious 
effects of this limitation and the utility of the discrete 
Fourier transform for transient signal analyses may therefore 
be diminished. Various data windows have been used in the 
past in an attempt to improve performance, but all involve 
tradeoffs in resolution and sidelobe levels that may 
unacceptably affect the analysis. [Ref. l:pp. 63-82] 

Recent appearance in the literature of several articles 
describing wavelet transforms and the related multiresolution 
analyses have created interest in determining their 
suitability for signal analysis applications; specifically, to 
overcome the time localization limitation of the Fourier 
transform. This thesis examines the potential of the wavelet 


transform for signal analysis purposes. 


II. REVIEW OF WAVELET THEORY 


A. DISCRETE WAVELET TRANSFORM 


Wavelets are families of functions, 


a p(X) =|al/?2 (=> =| ae (25 


generated by the dilations and translations of a single 
function. For digital signal processing applications it is 
necessary to discretize the parameter values and fix the 


initial dilation step a, and the translation step b, such that 


57 be Da 


Equation 2.1 becomes, 


h,, ,(x) =ag’*h(agx-nb,) m,neZ (ues 


with a=a,”" and b=nb,ag” 


The translation parameter is therefore dependent on the 
Gilation parameter. A large, positive m will contract the 
function hyo and cause the translation step to shorten while 
a large negative m dilates the function hyo and lengthens the 
translation step size. The inverse relationship ensures 


coverage of the entire range. [Ref. 2: pp. 909-910] 


Associated with the discrete wavelet is the discrete 
wavelet transform, T, which maps the function f to a sequence 
indexed by z°, a are 
(Try Ghee =a) [ Blagx-nb,) £ (x) dx ees 


ita 


flelt1F (hy) Rae < &, (2.4) 


and h has sufficient decay, and if T has a bounded inverse on 
meomebange, the set <h,> forms a frame for all square- 


integrable, one-dimensional functions f, 


£ (x) €L? (R) (235) 


The significance of a frame is that numerically stable 
algorithms exist which allow f to be reconstructed from the 
wavelet coefficients <h,,,f>. [Ref. 2: p. 911] 

Selection of appropriate h,ao, and bo is influenced by 
various factors. The desire to take advantage of previously 
published work and to minimize redundancy in the wavelet 
representation resulted in choices of ap = 2, bo = 1, and 


several h (discussed below) such that the h, constitute an 


orthonormal basis of compact support. 


B. MULTIRESOLUTION ANALYSIS 


Multiresolution analysis, as the name implies, describes 
the L* function f as a series of approximations of the 
function at different levels of resolution and consists of a 


family of embedded closed subspaces 


VC L*(R): |. SEV Reve eVelneVeraa (226) 
such that 
Vv, = {0}, Uv, = £7(8) ao 
and 
PU) 6 Va eae C2035) 


Also, there is a scaling function (x) €V, such that the 


dn (where ,,,(x) =27/(2"x-n)) constitute a basis for V, 
Vea span (@,,,) (2.9) 


Let Pf be the orthogonal projection of f onto V.. Pronmene 


preceding, it can be seen that lim_.P.f = f, for all feEL*(R). 


(Ref. 3:pp. 967-968] 


Considering our parameter choices and the decision to use 


orthonormal bases, if we define 


rea (en) C25 LO) 


mn 


then 


Peto: CaN) Ome 6 -V,, G22 71a.) 


Now define Qf =P f4f- Pf and WwW, as the orthogonal complement 


Grey, (W, 2 V,) so that V,,, = V,0W,,- Thee es Ss |) ehe 


orthogonal projection of any f € L*(R) onto WwW, the 


orthogonal complement of V, in V Pee eercters Caled 


m+ 


versions of WwW, , 


ig (9) SIC ng alee e232) 


and the W_ are orthogonal spaces which sum to L*(R) 


L?(R) = @ w,, Cr 3)) 


Similar to the V,, there exists in W, a wavelet function y 


such that 


Ae span [Wry] (2.14) 


where W,,(x) = 2%?y(2"™x-n) . If we define 
d_n ( Lae (2, 155 
then 
Of = )id (ea, “eave (2.16) 


Naturally, the function may be reconstructed from its 
projections onto the bases vector spaces. [Ref. 2:pp. 916-918] 
1. Properties of the Wavelet-Scaling Function Pair 


Following is a brief summary of the relationship 


between the @and ww families. Since 9,(x)E€V,cV, , it can 


be written as 


oe 


bo (x) = YP (by (x) 1b, (x-3)) 6, (x-Z) (ZF 


N=-e 


Let 


h(n) = [40 (%) 9 (2x-n) dx 
pnt (2.18 
= 92) 200 (a0) 1, (x-2)) / 
so that 
(x) = v2 YY Ain) >, (x-2) 
00 (2.19) 
= 2)> h(n) $,(2x-n) 
The Fourier transform of this equation is 
®(2w) = H(w) P(w) (2.20) 
where we have defined 
H(o) = )> h(n)e 7 (Ayes 


n=~0 


To satisfy the above relationships, the following properties 


must hold: 


IH#(O)| = 1 
(2-22) 
IW(w) |? + J¥(w+n) |? = 1 


From Equation 2.20 it can be seen that the @ can thus be 


derived if H(w) is known since 


®(w) = I H(2-Pw) (2.23) 
p=1 


Knowledge of @ can subsequently be used to find wy . Since 


Wo (x) € Wc V, , it can be written as 
Wino) = YE Wo (2) by (HF)? 4 (°F) G22) 
Let 


{Q 
S 
| 


= [Wo (x) by (2x-n) 
= (2.25) 


1 
> 2 


(wy, (x) 0, (x-=)) 


and by a similar approach to that used previously we can find 


W(x) = v2)7 g(n)o, Gea) 
i (2.26) 


2 » g(n) , (2x-n) 


n=-0 


The Fourier transform of Equation 2.26 is 


Y(2w) = G(w) ®@(w) C257) 


The following properties result: 


IG(0)| = 0 
IG(@)|? + IG(w+n)|* = 1 (2-28) 
H(w)G(w) + H(wt+t) G(wtnt) = 0 


Equations 2.22 and 2.28 are recognized as characteristic of 


"Conjugate" filters and are necessary to ensure orthogonality 


among the @ and wW families of functions. Equation 2.28 


ensures that each function of the @q (commonly referred to 


as the scaling function) is orthogonal to each function of the y 


family. This last condition results from the fact that 


V, 


m 


1 W 


m e 


[Ref. 4: pp. 16-23] 


An example of a function that fulfills the stated 


conditions is 


G(w) = eJ*Hl@tn) (2-219) 


from which we can find 


g(n) = (-1)'*A(i=n) (2-365 


Equation 2.30 is the defining relationship for describing 
wavelet-scaling function dependency in the application 


algorithms. [Ref. 3: p. 6799 


2. Decomposition and Reconstruction Algorithms 


Recalling that c,(n) = (p,,,f) , we can derive the 


decomposition recursion algorithn: 


Cm-1 a) Onnaae: f) 


= (p__, (x-27(""1) n) , £) 
=u, (x- 2") n) F(x) ax 


= fwvzy h(k) , (x-27 1) n-2-™k) ] £ (x) dx 
kK 


(2,315) 
= V2) ACK) [q(x-2-™ (2ntk) ) £ (x) dx 
k 
Shy 2) Dig) emanate) 
k 
Likewise we find 
Gn (1) = ae o (K)iceaia) (2.32) 
k 
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To determine the reconstruction recursion algorithm recall 


Pit = Pyylt+On yt 


= C27,.33)) 
1G ao ae 7 b= Cc =i (K) ® cmeayx* Du d “1 (k) W (m-1)k C2 .34) 

and substitute terms to get 
Cn = ans Pnf) (2.35) 


= p> Cc -1 (k) nn Pemar* Ds d -1 (k) mn Wee 


Using a change of variables the following can be found 


(ans Omi) = ¥2h(n-2k) (2.36) 
Ree - ¥2g9(n-2k) (2.37) 


so that 


Cnn) = V2[D> Cys (Kk) h(n-2k) +) d,., (kK) g(n-2k)] (2.38) 
k k 


The equations of this section demonstrate the pyramidal 
structure of multiresolution analysis and readily lend 
themselves to digital signal processing applications. [Ref 4: 
mm. 25-28 | 

The important conclusion is that, assuming knowledge 
of the g and h vectors, the c and d coefficients at any level 
can be completely determined from the c coefficients at the 
next higher level; also, the c coefficients at any level can 


1d! 


be determined from the c and d coefficients at the adjacent 


lower level. Notice that calculations within this pyramidal 


algorithm do not require explicit use of the @ andvy 


functions since interlevel relationships are defined entirely 
by the g and h vectors. Obviously, the same results would be 
reached if dilations and translations of the orthogonal basis 
functions were used directly at each level in the 
decomposition and reconstruction, but would be found at great 
computational cost because of the requirement to calculate 


inner products. 
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IIIT. COMPUTER IMPLEMENTATION OF MULTIRESOLUTION ANALYSIS 


A. INTRODUCTION 

Our goal is to investigate the use of multiresolution 
analyses to extract information from an input signal. In this 
case, a one-dimensional data sequence will represent the 
iA put . To avoid having to numerically evaluate any inner 
products we will equate the c,(n) coefficients with the data 


sequence, thereby constructing an auxiliary function f, with 
f=Lic,(n)o,(n), which clearly resides in V,. Since c,(n) 


entirely describes the input, it represents the highest 
resolution level possible and the apex of the pyramidal 
algorithm. The multiresolution framework previously described 
can now be used to decompose f (thus the data sequence Cy (N) ) 
into lower resolution, i.e. m < 0, approximation coefficients 
c,(n), and detail coefficients d,(n) . The reconstruction 
recursive equations can be used to recover the original f (and 
therefore the data sequence). A different multiresolution 
analysis exists for each scaling function/wavelet set. MThis 
thesis is limited to the Haar wavelet and to Daubechies' group 
of compactly supported wavelets. The h vectors for these 
wavelets were published in Reference 1 and are listed in 
Appendix H for convenience. Describing the h vectors is the 
common method of defining scaling function/wavelet sets 
because all other desired information can subsequently be 


derived from them. 


AB) 


Several computer programs were written to implement 
various multiresolution analyses. The programs were written 
in MATLAB™ to take advantage of the transportability of m- 
files between various operating systems. The capabilities of 
these programs, listed in Appendices B-H, will be described in 
general. Specific information can be obtained by referring to 
program comments in the appropriate appendix. 

The normal entry point into the collection of programs is 
the calling program wavelet.m (Appendix 8B). A brief 
description of the various possible options is presented upon 
entry and the user can make a decision based on the choices 
provided. Specifically, the user may decide to use the Haar 
wavelet, his own wavelet, or one from the Daubechies group of 
nine compactly supported wavelets and analyze with either a 
single phase or multiphase approach. 

The single phase approach implies the decomposition start 
point is the data start point, l.e. a "snapshot" of the entire 
data sequence is processed. The term multiphase refers to 
starting a decomposition at every possible data 
sequence/wavelet phase relationship, i.e. a "sliding window" 
approach, to maximize the signal analysis capability. As 
desired, wavelets are phase sensitive but feature extraction 
may be hindered if the "best" input phase relationship is not 
used. The multiphase analysis therefore contains multiple 
sets of coefficients, with each set containing all the 


information necessary for reconstruction. 
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B. SINGLE PHASE MULTIRESOLUTION ANALYSIS 

The second program, haarwave.m (Appendix C), does the 
"snapshot" analysis of an input data sequence using the Haar 
basis. Because of the simple nature of the Haar, it is easy 
to generate the actual approximation and detail function at 
each resolution level. The input data is viewed as a 
piecewise constant function (equating to the sample and hold 
operation) to allow easy inner’ product '= calculation. 
Representations of the properly weighted scaling function and 
wavelet at each level are provided to clearly indicate the 
dyadic relationship between levels. Additionally, the 
reconstruction of the original sequence can be viewed as well 
as a plot of the reconstruction error. Finally, distribution 
of the energy of each approximation and detail coefficient and 
the total energy of each resolution level are displayed. 
These energy distributions are the foundation for signal 
analysis. 

The third program, daubwave.m (Appendix D), allows the 
user to pick one of the wavelets from Daubechies group or to 
enter a valid user-defined h vector to define the basis set. 
Basically, the same output plots that were generated for 
wvhaar.m can be seen. Because of the complexity and 
irregularity of the basis functions, however, only the value 
of the coefficients which are generated for specific points 


will be plotted instead of the inner product representation. 


dbs, 


C. MULTIPHASE WAVELET ANALYSIS 

The last analytical program, multiphs.m (Appendix EE), 
allows the use of any of the previously defined wavelets with 
the multiphase algorithms. The different coefficient sets are 
displayed simultaneously, which is possible because their 
coefficient indices interleave without interference. Clearly, 
redundant information is provided in the sense that there are 
more coefficients than would be necessary for reconstruction 
of the data sequence over its interval of support. However, 
because of the different zero padding necessary for each set 
of coefficients a different Pf in V, is defined for each set. 
The user will therefore be able to see the P,f with the most 


distinctive set of coefficients to ease analysis. 


D. SUPPORTING PROGRAMS 

The program basisplt.m (Appendix  F) Supports’ the 
daubwave.m and multiphs.m programs by generating iterative 
plots of the scaling function and wavelet bases. It can also 
be called up directly and will use the vector "hcoeff" as the 
h coefficients to determine the @@ and wW based on a 
graphical recursion method. 

The functions wvinput.m (Appendix G) and daubdata.m 
(Appendix H) are called as necessary by any of the 
multiresolution analysis programs to provide selected input 
data signals or the h vectors from Daubechies' group, 


respectively. The input signal options consist of a sinusoid, 
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a pseudo-noise (PN) sequence, a sinusoid modulated with a PN 
sequence, and any of the above corrupted by noise. Various 


parameters can be modified to satisfy the user's needs. 


AL 


IV. SIGNAL ANALYSIS 


All plots associated with this section can be found in 


Appendix A. 


A. SINGLE PHASE ANALYSIS OF A SINUSOIDAL WAVEFORM 

The first signal to be analyzed will be the sinusoid shown 
in Figure 1, along with its level 0 approximation. Figures 2- 
5 show the Haar multiresolution analysis, stopping at level -4 
since there is no energy at any lower level for this 
decomposition. The reconstruction plots (not shown) are 
indistinguishable from the decomposition plots because the 
maximum reconstruction error is 15 orders of magnitude below 
signal levels. The energy contained at each resolution level 
for both approximation and detail coefficients can be seen in 
Figure 6. Clearly, the maximum energy change occurs in the 
detail coefficients at level -4 and the energy in the 
approximation coefficients at every level is the sum of the 
energy of the detail and approximation coefficients from the 
level below, as expected. Figures 7 and 8 provide a clear 
summary of exactly where the signal and filter best match with 
respect to sample number and resolution level. In this 
example, the Haar decomposition highlights the phase changes 
associated with the switching between positive and negative 


half-cycles because of the signal/wavelet phase coherence. 
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The decomposition is repeated with a different 
multiresolution analysis based on the 6-coefficient (third 
order) Daubechies wavelet. Only the coefficient energy plots 
are shown in Figures 9-11 for comparison with the Haar plots. 
It can be shown that the time index of coefficients goes 
beyond the original support length of the data (potentially 
out to sample 257 in this example, although the plots were 
truncated at sample 90 for display purposes and because very 
little energy 1s associated with samples beyond 90). Although 
the coefficients outside the signal support range are 
generally very small, they are necessary for completeness. 
Computing these coefficients requires affixing zeros to the 
sequence, most significantly at the trailing edge. These 
"edge effects" are minimized at the leading edge (no 
coefficients are generated for negative time indices) by 
shifting the h vector so that nonzero values may initially 
exist only over [-(m-2),...,0,1] and then assigning the 
coefficient value to the index of the second term from the 
right (initially 0). This is intuitively satisfying from a 
causality viewpoint and also aids interpretation by only 
having coefficients outside the data sequence on one side. 
Leading edge effects consist of the influence of leading zeros 
necessary for filter coefficients with negative time indices. 
Note that edge effects increase as the resolution level 
decreases because lower resolution coefficients are influenced 


by a larger number of data points. The extent of edge effects 


Abe 


is also dictated by the number of h_ coefficients. 
Coefficients beyond the data sequence support will not have 
time indices greater than [(2 °—1l) (number of h coefficients - 
1)) for any of the decomposition programs and often 
significantly fewer (e.g., if the data sequence length is a 
power of 2 the Haar decomposition will not generate any terms 
beyond the data length, as can be seen in the previous plots). 

The next set of plots (Figs. 12-18) show the phase 
sensitivity of the decomposition when the input sinusoid phase 
is shifted from 0 to 45 degrees. The energy distribution is 
no longer as concentrated as in the previous set of plots and 
interpretation is consequently more difficult. Changing the 
sampling frequency relative to the sinusoid frequency will 


have a Similar effect. 


B. SINGLE PHASE ANALYSIS OF A BINARY PHASE SHIFT KEYED SIGNAL 

The Binary Phase Shift Keyed signal shown in Figure 19 
will serve to highlight a deficiency in the single phase 
approach. Notice with the Haar wavelet in Figures 20-22 that 
the decomposition is exactly the same as the pure sinusoid, 
i.e., because of the phase alignment we cannot discriminate 
between the two signals' coefficient energy plots. The 
Daubechies 6-coefficient wavelet decomposition (Figs. 23-25) 
does appear significantly different than the sinusoid 
decomposition but there is no clear indication of every phase 


reversal. Thus, the results are ambiguous. 
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C. MULTIPHASE ANALYSIS OF A BINARY PHASE SHIFT KEYED SIGNAL 

Since these responses would be unsuitable for signal 
analysis purposes the multiphase approach (described 
previously) was developed. It can be viewed as starting 
decompositions at each of the 2" phases of each level so that 
the most distinctive representation can be used for analysis 
regardless of initial signal phase. The multiphase plots are 
shown in Figures 26 and 27 for the Haar and in Figures 28 and 
29 for the Daubechies 6-coefficient filter. There is clear 
indication in both sets of plots of the signal's phase 


maversions. 


D. MULTIPHASE ANALYSIS OF A TRANSIENT SIGNAL 

Finally, the transient signal of Figure 30 was 
investigated. Use of the "zoom" feature in the programs was 
necessary to clearly see the response. Figures 31-34 show the 
responses associated with the Haar and Daubechies 6- 
coefficient wavelets that have been the standard throughout 
this thesis. Higher order Daubechies (more regular) wavelets 
were also experimented with in an attempt to best fit the 
Signal and it can be seen that the 18-coefficient (ninth 
order) wavelet used for Figures 35 and 36 does a good job of 
clearly indicating the presence of the signal and of 
maximizing the concentration of energy into only a few points 
in a single resolution level. This good "match" of the signal 
with the wavelet filter could be used for’ signal 
Classification purposes as well as detection if ana priori 


PAD, 


selection of the wavelet filter based on a similar known 


signal had occurred. 
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V. CONCLUSIONS 


Wavelet-based multiresolution analysis is another tool for 
the signal analyst and great potential exists for its use in 
various applications. Signals that can only be represented in 
the frequency domain with large numbers of significant terms 
provide special motivation for wavelet-based decomposition, as 
can be seen in the transient Signal example above. Signal 
detection and pattern recognition through the use of "matched" 
wavelets may also prove useful, especially in areas where 
analysis at different resolution levels (1i.e., possible signal 
Gilation/contraction) is required. The extent of the ultimate 
usefulness and popularity of wavelets will become known as 
others gain knowledge of basic wavelet characteristics and 


incorporate multiresolution analyses into their research. 


23 


APPENDIX A 


SIGNAL ANALYSIS PLOTS 
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Figure 8. Haar Detail Coefficient Distrieucren 


JOAY] UuOINOssy 

I 0 is Ga og p- ¢- 9- ie 
0 
v0 
9°0 


80 


(Gl 
SJUDIDIJJIOZ [!VJaq JO Ad1aUz] pozi[peulION 


[OAT] UONN[OSaYy 


S- 2° Le 


[- C oi a Bore a eee 0 
c0 
0 
90 
8°0 


(a 
SUdIDIJJIOD UOHLUNXOoIddy jo Ad1auq pazlpeulION 


Daubechies Resolution Level Energy Distribution 


Figure 9. 


33 


Figure 10. Daubechies Approximation Coefficient Distribution 


Figure 11, Daubechies Detail Coefficient Distribution 
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Figure 14. Haar Approximation Coefficient Distribution 


Figure 15. Haar Detail Coefficient Distribution 
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Figure 17. Daubechies Approximation Coefficient Distribution 


Figure 13. Daubechies Detail Coefficient Distribution 
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Figure 22. Haar Detail Coceti retenme Distribugven 
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Figure 24. Daubechies Approximation Coefficient Distributi#em 


Figure 25. Daubechies Detail Coefficient Distribution 
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APPENDIX B 


CALLING PROGRAM FOR ACCESSING MULTIRESOLUTION ANALYSES 


PROGRAMS AND FUNCTIONS 
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in the comments below. 


I 
— SS 
se 


specific needs. 


Program options include: 


a 
Haar wavelet 


a. 
ae 
3. Selection of Input Data 


1) Sine wave 


b. User~-input data vector 


™ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ™~ ~ ~ ~~ ow 


Calling program for wavelet multiresolution analysis. 
and options of the entire set of programs are presented, in general terms, 


the program (after prompts and questions, 


The capabilities 


to view next plot, etc.) 


Selection of Basis Functions 


b. Daubechies group of compactly supported wavelets 
c. User-input wavelet ("h" coefficients) 


Generation of Scaling Function/Wavelet plots 


a. User-defined data vector (with optional noise background): 


2) Pseudo~-noise (PN) sequence 
3) Sine wave modulated by PN sequence 


<Returm> for more <2. 


N2 


4. 


Decomposition algorithm based on: 


a. 


Data "Snapshot" - Fixed start point; 


one decomposition; 
Single phase approach 


Da 


Do you want to see: 


se 


ba ~ ™~ ™ ~ ~ ™~ ~ ™~ ~ ~ ™~ ~ ™ ~ ~ ~ ~. ~ = 


oe 


co I es 
disp(b) 
disp(Nl) 
pause 
cle 
disp(b) 
disp(N2) 
disp(b) 
desire 


input(’Answer 1, 


"Sliding" Data window - Floating start point; decomposition 


for each start point; multiphase 
approach; discrete approximation to 
continuous wavelet transform 


Haar wavelet data "snapshot" decomposition, 
(input data treated as a piecewise constant function, 
allowing inner products to be calculated and displayed) 
2. Daubechies or user wavelet data "snapshot" decomposition, or 


"Sliding" data window wavelet (any type) decomposition? 
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This program will perform multiresolution analysis on input data’ 
through the use of various orthonormal bases of compact support. 
It is interactive in nature and will allow you the opportunity 
to select among several different options to best support your 
Use the Return key to control movement through 


/ 
é 
é 
s 
/ 
/ 
s 
f 
rf 
f 
/ 
f 
f 
/ 
s 
s 
/ 
, 
rf 
, 
f 
f 


_—— 


~ »~ ~ ws ~~ ww ww w %*% ~~ = & % ~ ~ ~ ~~ Ss % 


=e 


“=e 


if desire == 1, 
haarwave 

elseif desire == 2, 
Gaubwave 

else 
multiphs 

end 
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APPENDIX C 


PROGRAM FOR SINGLE PHASE ANALYSIS USING HAAR WAVELET 
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% This program performs single-phase Haar wavelet decomposition on a 
$ user-supplied or program-generated input waveform. Representation is 
% made with bar graphs to show actual inner product calculation (input 
¢ data is taken as a piecewise constant function). The input data sequence 
$ will be zero-padded to accomodate the calculation of all possible non- 
% zero coefficients and approximation ("c") and detail ("d") coefficients 
$ may be computed for sample numbers beyond the data endpoint. The user 
% should consider these "edge effects" in his analysis. 
% Determine if user wants to input an external data vector or desires to 
* build one through the program. 
cic 
ots oe 
Q1 = [{’Do you wish to use: f 

‘1. your own MATLAB formatted row vector, or : 

‘ 2. a program generated vector from the following menu?’ 

: - Sine wave ‘ 

z - Pseudo-noise (PN) sequence 

‘ - Sine wave modulated by PN sequence cole 
disp(b) 
disp(Q1) 


Pick = input(’Answer 1 or 2: ‘); 


merck —— 1}, 
% Read in user’s input vector 
Ni = [’ Note: If the number of data points is a power of 2, the’ 
’ results are easier to interpret because there are not ’” 


‘ any “edge effects". alee 
disp(b) 
G@isp(N1) 
Wavedata input(‘/What is the name of your input vector? ‘'); 


sampfreq = input(’What was the sampling frequency (Hz)? '); 
Tsample = 1/sampfreq; 


else 
{[Wavedata,Tsample] = wvinput(Pick); 
end 
% Decomposition algorithm n(0)—055,- h(1)=0.5, 9g (0)=0.5,00{1)=-0.5 


datingth = length (Wavedata) ; 
numrows = ceil(log(datlngth) /log(2)); % find number of resolution levels 
numpts = 2“numrows; Newdata = zeros(1,numpts); 
Newdata(l:datlngth) = Wavedata; 
dq = zeros(numrows,numpts); ¢ = zeros(numrows+1,numpts) ;} 
detail = zeros(numrows,numpts); approx = zeros(numrows,numpts) ; 
c(numrows+1,:) = Newdata; 
top = numpts; 
const = 1/sqrt(2); % normalization constant * coefficient magnitude 
ctr = numrows; 
Give ctr >= l, 
n = numrows-ctr+l; 
shift = 2“n; 
shiftl = shift/2; 
top = ceil(top/2); %* number of points at this level 
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for k = l:top % get "c" and "dl" coctticrenes 
jump =shift*k; 
indexl = jump-(shift-1); 
index2 = jump-(shifti-1); 
firsttrm = (¢(etr+l, 1 ndex1)-- 
if index2 <= numpts, 
secndatrm = c(ctrt+l,index2) ; 
else 
secndtrm = 0; % zero padding 
end 
a(ctr, index1) firsttrm - secndtrm; 
C(Ctr, index?) firsttrm + secndtrm; 
for j) = 0:(shift-1) % build matrices for display purposes 
if 30 < Shier 
detail(ctr,indexl+j) = ad(ctr,indexl); 


not 


else 
detail(ctr,indexl+j) = -d(ctr,indexl); 
end 
approx (ctr, indexl+)]) = c(ctr,indexi): 
ena 

ena 
Gi(Ctr,:) = Const*a(ctrr.:):; 
G(Ctr,:) = const*c(etr,.<); 
normlize = sgqrt(2*shift); % normalization constant for this level 


Getail(ctr,:) 
approx (Ctr, <) 
CLE = €vr-1; 


Getail (ctr, :) /normlize- 
approx (ctr) /nerm)ize-, 


end 


6 Plot "c" and "a" coefficients by resolution level 


indxtod [0.5:numpts-0..5)? 

plotmin 1.2*min(Wavedata); plotmax = 1.2*max(Wavedata) ; 

if ((plotmin==0) &(plotmax==0)), plotmax = 0.5; end 

1f plotmin > 0; plotmin = 6+ end 

if plotmax < 0, plotmax = 0; end 

v=[0,numpts,plotmin,plotmax]; 

axis(v); 

bar (indxto0,c(numrows+1,:) ) 

title(’Approximation function to Input Signal at resolution level 0’) 
xlabel([’Sample number "n" (Time = n * ’,num2str(Tsample),’ sec) ’)) 
pause 

eld 


outptdet = zeros(1,1); outptapp = zeros(1,1); 
for KK = numrows:=1:1 
level = K-numrows-1; 
step = (2*(-level))/2; 
for j} = l:numpts/step 
outptdet(j) = detail(k, (j~-1)*step+1) ; 
end 
stepa = 2*step; 
indxto0a = [stepa/2:stepa:numpts-stepa/2]; 
for } = l:numpts/stepa 
outptapp(j) = approx(k, (j-1)*stepa+l) ; 
end 
axlis(v); 
Ve eee, 
subplot(211),bar(indxto0a, outptapp) 
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else 
indxto0a = [O:numpts]); 
approxl = approx(1,:);approxl(numptstl) = approx(1,numpts) ; 
subplot (211),plot(indxto0a, approxl1) 
end 
title({’Approximation function at resolution level ’,num2str(level) }) 
xlabel({’Sample number "n" (Number of samples = ’,num2str(numpts),’)’)) 
axis(v); 
subplot(212),bar(indxto0, outptdet) 
title([’Detail function at resolution level ’,num2str(level)}) 
xlabel([{’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’}) 
pause 
Cig 
clear outptapp outptdet indxtoo 
indxtoo = indxto0a; 
clear indxto0da 


subplot 


% Reconstruction algorithm 


Elc 
recmpout = zeros(1,1); 
disp(b) 
N2 = [’ Level-by-level Recomposition can be observed if desired. 
’ The Recomposition algorithm starts with the lowest level ’ 
’ Approximation Function and successively adds in the Detail’ 
’ Function to obtain the next higher level Approximation. ees 
disp(N2) 
disp(b) 
ckalgthm = input(’Do you want to see the Recomposition (Y or N)? ’,’s’); 
Mmerckalgthm == ‘’Y’, 
level = -numrows; 
axis(v); 


plot (indxto0, approxi) 
title({’Level ’,num2str(level),’ Recomposition’ )) 
xlabel(’Sample number’) 
pause 
Clg 
Pecomp = approx(1,:); 
for k = 1:numrows 
level = level+1; 
recomp = recomp+detail(k,:); 
step = 2%(-level); 
indxtooO = [step/2:step:numpts-step/2}); 
for j = l:numpts/step; 
recmpout(j) = recomp((j-1)*steptl); 
end 
axis(v); 
bar (indxto0o, recmpout) 
title([’Level ’,num2str(level),’ Recomposition’ }) 
Xlabel(’Sample number’) 
pause 
cig 
end 
bar (indxto0, Newdata-recmpout) 
title(’Recomposition Error’) 
xlabel(’Sample number’) 
pause 
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Coefficient energies are used as a measure of response 


C=C. 727 

C= 70.727, 

Enrgytot = sum(c(numrows+1,:)); 

Enrgycro = zeros(1,numrows+1); Enrgydro = zeros(1,numrows- 1) 

for j} = l:numrows % find energy in each resolution level 


yy / ENEGY Con, 
)) JEnECy coc. 


Enrgycro(j) = (sum(c(j,:) 
Enrgydro(j) (sum(d(j,:) 

end 

if Enrgytot == 0, 
Enrgycro (numrows+t1) 

else 
Enrgycro(numrowst+1) = 1; 

end 

Enrgydro(numrows+1) = QO; 

lvl = [=(numrows): 1:0]; 

Vv = [=(numrows+t?) 20.1.2); 

axis(v); 

subplot(211),bar(llvl, Enirgyere) 

title(’Normalized Energy of Approximation Function vs. Resolution Level’) 

xlabel(’Resolution Level’) 

axis(v); 

Subplet(212),bar(ivi,Enreydro) 

title(’Normalized Energy of Detail Function vs. Resolution level’) 

xlabel(’Resolution Level’) 


0; 


co oe ee oe ee ee ee ee es es ee ee ee es ee es ee ee ee ee i ie ee ie ie ie eee eee ee eee ee ee ee ee ee ee ee eee ee eee eS eee ese 


Display mesh and contour plots of coefficient energy with optional "zoom" 
% capability to aid analysis 


strtsamp = 1; endsamp = numpts; strtrow = 1; endrow = numrows; zoom = l? 


highres = -1; lowres = -numrows; 
while zoom == 1, 
rangea = [{-highres-1:-lowres]; ranged = [-highres:-lowres]; 


indxtooO = [strtsamp-1:endsamp-1]; 
mesh(c(strtrow:endrow+l1,strtsamp:endsamp) ) 

title(’Energy in Approximation Coefficients using Haar basis’) 

pause 

contour(c(strtrow:endrowt+1,strtsamp: endsamp) ,10,indxto0,rangea) 
title(’Contour Map of Approximation Coefficient Energy Distribution’ ) 


xlabel([’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number (= -Resolution Level)’) 
pause 


mesh(d(strtrow:endrow,strtsamp: endsamp) ) 

title(’Energy in Difference Coefficients using Haar basis’) 

pause 
contour(d(strtrow:endrow,strtsamp:endsamp) ,10, indxto0, ranged) 
title(’Contour Map of Difference Coefficient Energy Distribution’ ) 


xlabel(([‘Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number (= -Resolution Level) ’) 
pause 
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ele 
Siowmoce.= i1nput« Would you Pike to "zoom" in on a section (Y or N)? ‘,’s'); 
if showmore == ‘Y’, 


aisp(b) 
Gisp(’ You will set the display sample number and resolution level limits. ’) 
adisp(b) 
aisp([” Sample range O:’,num2str(numpts-1) }) 
Gispet! Resolution range —~1:’,num2str(~numrows) }) 
revu = input(’Need to see the original contour plots again (Y or N)? ’,’s‘); 
if revu == 'Y’, 
Goncour(c,10,[O0:numpts—1],[0:numrows ]) 
title(’Contour Map of Approximation Coefficient Energy Distribution’) 


xlabel([{’Sample number "n" (Time = n * ‘,num2str(Tsample),’ sec)’)) 
ylabel(‘’Decomposition Number (= -Resolution Level) ’) 
pause 


eencour(a,210,([0snumpts—1)], (1: numrows]}) 
title(‘’Contour Map of Difference Coefficient Energy Distribution’) 


Xlabel([‘Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number (= -Resolution Level) ’) 
pause 

end 


strtsamp = input(’Sample start point? ‘’)+1; 

endsamp = input(’Sample endpoint? ’)+1; 

highres = input(’Highest resolution level (least negative)? ’); 
lowres = input(‘’Lowest resolution level (most negative)? ’); 
endrow = numrowstit+highres; 

strtrow = numrows+l+lowres; 


else 


zoom = 0; 


end 
elg 


end 
elec 
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APPENDIX D 


PROGRAM FOR SINGLE PHASE ANALYSIS USING 


DAUBECHIES OR USER WAVELET 
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This program performs single phase decomposition of the input waveform 
with either a Daubechies wavelet of user-defined order or with a wavelet 
provided by the user. The input data sequence will be zero-padded to 
accomodate the calculation of all possible non-zero coefficients and 
approximation ("“c™") and detail ("d") coefficients may be computed for 
sample numbers beyond the data endpoint. The user should consider these 
"edge effects" in his analysis. 


O° &© of oO cl? ol oO 


DD i et te ini 


% Input the "h" coefficients 


ele 

ee LY): 

disp(b) 

Dil = {’ Would you like to use: i 
j 1. your own decomposition coefficients (use wvhaar.m if you desire’ 
ij to use the Haar basis set)’ 
2. those associated with Daubechies compactly supported wavelets? ‘’]; 

disp(D1) 

Morek = input(’Answer 1 or 2: ‘); 

disp(b) 

if hpick == 1, 


disp(’ Input as a row vector with sum of coefficients normalized to "1".7’) 
hcoeff = input(’What is the name of your decomposition coefficient vector? ’); 
choice = 1; 

else 
choice = input(’What is the desired wavelet order (2-10 are available)? ’); 


{heoeff) = daubdata(choice) ; 
end 


hlength = length(hcoeff) ; 


if hlength == 0 
disp(b) 
disp(’ERROR: Wavelet order selected is outside allowable range. ’) 
break 


% Determine if user wants to see plots of basis functions 


ckplt = input(’Do you want to see the Scaling Function/Wavelet (Y or N)? ’,’s’); 
if ckplt == ‘y’ 

basisplt 
end 


* Get the input data vector 


ee 
disp(b) 
Q1 = {’Do you wish to use: , 


‘1. your own MATLAB formatted row vector, or 4 
2. a program generated vector from the following menu?’ 
- Sine wave ; 
- Pseudo-noise (PN) sequence 
- Sine wave modulated by PN sequence ies 


f 


~ ~ % %& 


V1 


disp(Q1) 
Pick = input(’Answer 1 or 2: ’); 


1f Pick == 1, 
% Read in user’s input vector 
Wavedata = input(’What is the name of your input vector? ’); 
sampfreq = input(’What was the sampling frequency (Hz)? ’); 
Tsample = 1/sampfreq; 


else 
[Wavedata,Tsample] = wvinput(Pick); 
end 
Ge eee ewan an eee es escee ce eee a acnnase eee == ee 
% "g" coefficients are derived from the "h" coefficients 
hcoeff = sgrt(2)*hcoeff; 
gcoeff = fliplr(hcoeff); 
for j = 2:2:hlength 
gcoeff(j) = -gcoeff(j); 
end 
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% Decomposition algorithm 


datingth = length (Wavedata) ; 
numrows = ceil(log(datlngth)/log(2)); % find number of resolution levels 


Lastpts = zeros(1,numrows+1); Shift = ones(1,numrows+1) ; 
Lastpts(numrows+1) = datingth; lastpt = datilngth; 
for k = numrows:-1:1, % find number of points required at each level 
evnorodd = rem(lastpt,2); 
if evnorodd == 0, 
lastpt = lastpt/2+ceil((hlength-2)/2); 
else 
lastpt = (lastpt-1)/2+fix(hlength/2); 
end 
Lastpts(k) = lastpt; 
Shift(k) = 2% (numrows-k+1) ; 
end 


hhalf = ceil(hlength/2); 
adprime = zeros(numrows, Lastpts (numrows) thhalf) ; 
cprime = zeros(numrows, Lastpts (numrows)+hhalf) ; 


clastvect = zeros(1,datingth+2*hlength=3) ; 
clastvct(hlength-1:datlngththlength-2) = Wavedata; 
ctr = numrows; 
while. Ctr =>. 0; 

lastpt = Lastpts(ctr); 

nwevctr = zeros(1i,lastpt); 

nwdvctr = zeros(1,lastpt); 


for k.= lJ leastae % coefficients calculated, by resolution level, 
startpt = 2*k-1; % using convolution operation with shifts of 2 
endpt = 2*k+thlength-2; % (downsampling) 


NWCVCtr(k) = heoefE*clastVet(startpe:endoey, 


nwdavectr(k) gcoeff*clastvct(startpt:endpt) ’; 
end 
clastvct = [(zeros(1,hlength-2) nwevctr zeros(1i,hlength-1) ]); 
eprime(ctr, i: lastpt) = nwevcer: 


V2 


dprime(ctr,l:lastpt) = nwdvctr; 
cere CEY=<1> 
end 


% Coefficient energies are used as a measure of response 


cenergy = cprime.%2; 
denergy = dprime.%*2; 


The next plot will show the normalized energy distributions of the’ 
approximation and detail coefficients as a function of resolution 
level. Based on this information, you will select the minimum 
resolution level to be displayed on all future plots. 


Nl = [ 


"Edge effects" can cause coefficients to be generated for sample 
numbers beyond the endpoints of your data sequence (the lower the 
resolution level, the greater the sample number required; therefore, 
zero padding of the original data sequence is required). 


‘ 
f 
é 
f 
é 
‘ 
é 
é 
é 
This program translates the indices of the "h" and "g" vectors to ’ 

{[ -(length of vector-2),1 ]. As a result, coefficients will exist , 
for sample numbers beyond the last data point (n=N), but will not ; 
exist for sample numbers prior to the first data point (n=0). 
These coefficients are necessary for completeness, but may be Z 
difficult to interpret because the number of original data points : 
used in their calculation can be a small percentage of the total : 
considered (the rest are "0"). Similarly, coefficients calculated ‘ 
for sample numbers near the beginning of the sequence will also be Z 
affected by leading zeros affixed to the data sequence. 
c 

‘ 

é 

é 


Limiting the sample numbers required for display may assist inter- 
pretation; select the lowest resolution level containing significant 
energy for this effect. <Return> 


~~ * » & ® © ~ ~& * * ®S ~~ %S ~~ BS B® SS SX BH BY BY BW ® B® 


Ihe 


elec 
disp(b) 
disp(N1l) 
pause 


originlE = zeros(1,datlngth); originlE = Wavedata.“%2; 

Emagytot = sum(originlE) ; 

Enrgycro = zeros(1,numrows+1); Enrgydro = zeros(1,numrows+1) ; 

ckerr = zeros(1,numrows) ; 

for k = 1:numrows % find energy in each resolution level 
Enrgycro(k) = sum(cenergy(k,:)); 


Enrcydro(k) sum(denergy(k,:)); 
end 
Pepenrgytot == 0, 
Enrgycro(numrows+1) = 0; 
else 
Enrgycro(numrows+1) = Enrgytot; 
eee =—ENBovenrO,ENragycOt; Enrgydro = Enrgydro/Enrgytot; % normalize 
en 
Enrgydro(numrowst+1l) = 0; 


ckenergy = Enrgycro(il:numrows)+Enrgydro(1:numrows) ; 

* Energy balance check. The total energy at any level should equal the 
* energy in the "c" coefficients in the next higher level. 

for kK = l:numrows 


13 


ckerr(k) = abs(Enrgycro(k+1)-ckenergy(k)); 

if ckenergy(k) > 10*(-6)*Enrgytot, 

errenrgy = ckerr(k)/max(Enrgycro(k+1),ckenergy(k) ); 

if errenrgy > 10~(-6), 
Ele 
disp(b) 
disp(’ ERROR: Energy check between levels does not balance. ’) 
Gdiep( If you entered your own "h" vector, recheck its validityoas 
break 

end 

end 

end 


ivi = (=Qiumrows):120)]. 

V = (-(numrows+1) 71071221, 

axis(v); 

Subplot(211), bar (ivi; enrgyero) 

title(’Normalized Energy of Approximation Coefficients’) 
Xxlabel(’Resolution Level’) 

axis(v); 

subplot (212), bar (ivi, Enrgydro) 

title(’Normalized Energy of Detail Coefficients’) 
xlabel(’Resolution Level’) 

subplot 

pause 


Output the number of points required to properly display coefficients at 
each resolution level (determined by "edge effects") 


Cle 
disp(b) 
N2 = [’ Listed below are the number of samples required to properly display’ 
‘ each Resolution Level. “le 
disp(N2) 
disp(b) 


disp(’ Resolution Level Number of samples’) 
for k = l:numrows 


vetrindadx = numrows-k+1; 
numsamps = (Lastpts(vctrindx)-1) *Shift(vctrindx) +1; 
@isp Gl noumnesterc- i) ‘, num2str(numsamps) )) 
end 
nmbrlivl = -input(’What is the lowest Resolution Level you desire to see? ’); 


* Plot “c" and "d” coefficients by resolution devel 


indxtoo = (O0:datingth-1); 

plotmin = 1.2*min(Wavedata); plotmax = 1.2*max(Wavedata) ; 
if ((plotmin==0) &(plotmax==0)), plotmax = 0.5; end 

if plotmin > 0, plotmin =s0; end 

if plotmax < 0, plotmax = 0; end 

v = ([0,datlingth-1,plotmin, plotmax); 

axis(v); 

plot(indxto0,Wavedata, ’*’) 

title(’Approximation Coefficients at resolution 0”) 
xlabel({’Sample number "n" (Time = n * ’,num2str(Tsample),’ sec) ’}) 
pause 


lastrow = numrows-nmbrlvl+l1; 
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clmnsize = (Lastpts(lastrow) -1)*Shift(lastrow)+1; 
ee——szerostnmbrlVvi+1,clmnsiZe) ; 
G(mmorivi+i,1:datingth) = originlE;: 
ad = zeros(nmbrlvl,clmnsize) ; 
ctr = numrows; ctrl = nmbrivl; 
while ctr >= lastrow, 
sampval = zeros(1,Lastpts(ctr)); 
for k = 1:Lastpts(ctr) 
index = Shift(ctr) *(k-1)+1; 
c(ctrl,index) = cenergy(ctr,k); 
d(ctri,index) = denergy(ctr,k); 
sampval(k) = index-1; 
end 


ploumin = 1.2*min(min(cprime(ctr,:)) ,min(dprime(ctr,:))) 
plotmax = 1.2*max(max(cprime(ctr,:)),max(adprime(ctr,:))) 
if ((plotmin==0) &(plotmax==0)), plotmax = 0.5; end 
ieweloumin > 0, plotmin = @; end 

if plotmax < 0, plotmax = 0; end 


v = [0,max(sampval) ,plotmin,plotmax); 
n = numrows-ctr+1; 
axis(v); 


Euoploc( 211), plot (sampval,¢cprime(ctr,1:Lastpts(ctr)),’*’) 
title({’Approximation Coefficients at resolution level ’,num2str(-n)j) 
xlabel([{’Sample number "n" (Last sample = ’,num2str(max(sampval)),’)’)) 
Sumpleot(2i2) ,pleot(sampval,dprime(ctr,1:Lastpts(ctr)),’*’) 

title(({‘Detall Coefficients at resolution level ’,num2str(-n)}) 


Xlabel([{’Sample number "n" (Time = n * ’,num2str(Tsample),’ sec)’)) 
pause 
elg 
Pee — CLr—-1l; ctrl = ctri-1; 
end 
subplot 


% Reconstruction algorithm 


N3 = [’ Recomposition of the original data sequence can be checked if , 
’ desired. The recomposition algorithm starts with the approximation’ 
.' coefficients at the lowest resolution level previously selected ‘ 
’ and successively "adds" in the detail coefficients of each level 
’ until the original resolution of the input sequence is obtained. a 
elc 
Gisp(b) 
Gisp(N3) 


ckalgthm = input(’Do you want to see the Recomposition (Y or N)? ’,’s’); 
teeckalgthm == ‘’Y’, 
hrecomp = fliplr(hcoeff); grecomp = fliplr(gcoeff); % invert in time 
cstart = cprime(lastrow,:); 
ctr = lastrow; 
while ctr <= numrows; 
lastpt = Lastpts(ctr); 
Ckvectr = zeros(2,lastpt); 
for k = 1:fix(hlength/2) % compute higher level "c" coefficients 
ckvctr(1,:) = ckvctr(1,:)+hrecomp(2*k) *cstart(k:lastpt+k-1)+.. 
grecomp (2*k) *dprime(ctr,k:lastpt+k-1); 
SivGen(2, .)je—sCkVCEr(2,:)+hrecemp(2*k-1) *estart (k: lastpt+k-1)+.. 
grecomp(2*k-1)*dprime(ctr,k:lastpt+k-1) ; 
end 
if rem(hlength,2) -= 0, 
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ckvectr(2,:)=ckvctr(2,:)+hrecomp(hlength) *cstart(hhali: lastptahnalt—=ieeee 
grecomp(hlength) *dprime(hhalf:lastpt+hhalf-1) ; 
end 
ckvctr = reshape(ckvctr,1,2*Lastpts(ctr)); 
CS€art = ckvctr; 
lastsamp = (Lastpts(ctr+1)-1)*Shift(ctrtl) ; 
indxtooO = (0:Shift(ctr+1):lastsamp]; 
plotmin = 1.2*min(ckvctr); plotmax = 1.2*max(eckvetrEe 
if ((plotmin==0) &(plotmax==0)), plotmax = 0.5; end 
if plotmin > 0, plotmin = 0; end 
if plotmax < 0, plotmax = 0; end 
v = [(0,lastsamp, plotmin, plotmax]}; 
axis(v); 
plot(indxtoo0, ckvctr(1: Lastercs(ctr+i)) = 
title([{’Level ‘’,num2str(ctr-numrows) ,’ Recomposition’}) 
xlabel([{’Sample number "n" (Last sample = ’,num2str(lastsamp),’)’)) 
pause 
clg 
Ctr = CUrtl; 
end 
axis; 
plot (indxto0,ckvctr(1:datlngth) -Wavedata, ’*’) 
title(’Recomposition Error’) 
xlabel(([’Sample number "n" (Last sample = ’,num2str(lastsamp),’)’]) 
pause 


% Display mesh and contour plots of coefficient energy with optional "zoom" 
% capability to aid analysis 


strtsamp = 1; endsamp = clmnsize; strtrow = 1; endrow = nmbrlvl; zoom = 1; 


highres = -1; lowres = -nmbrlvl; 
while zoom == 1, 
rangea = [-highres-1:-lowres]; ranged = [-highres:-lowres]; 


indxtoO = [strtsamp-1l:endsamp-1]; 
mesh(c(strtrow:endrow+l1,strtsamp: endsamp) ) 
title(’Energy in Approximation Coefficients’) 
if choice ~= 1, 
xlabel([’using Daubechies Wavelet of order ’,num2str(choice) }) 
end 
pause 
contour(c(strtrow: endrow+1,strtsamp:endsamp) ,10,indxto0,rangea) 
title(’Contour Map of Approximation Coefficient Energy Distribution’) 


xlabel([{’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number (= -Resolution Level) ’) 
pause 


mesh(d(strtrow:endrow,strtsamp:endsamp) ) 
title(’Energy in Detail Coefficients’) 
if choice -= 1, 

xlabel([’using Daubechies Wavelet of order ’,num2str(choice) }) 
end 
pause 
contour(d(strtrow: endrow,strtsamp:endsamp) ,10, indxto0, ranged) 
title(’Contour Map of Detail Coefficient Energy Distribution’) 
xlabel({’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec) ’}) 
ylabel(’Decomposition Number (= -Resolution Level) ’) 
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pause 


ELC 
showmore = input(’Would you like to "zoom" in on a section (Y or N)? ’,‘’s‘'): 
if showmore == ‘Y’, 
disp(b) 
disp(’ You will set the display sample number and resolution level limits.’) 
disp(b) 
aispe| “ Sample range O:’,num2str(clmnsize-1) }) 
dasp([” Resolution range -1:’,num2str(-nmbrilvl) j) 
revu = input(’Need to see the original contour plots again (Y or N)? ’,’Ss’); 
if revu == ‘Y’, 


eont-our(c, 10, {[0:clmnsize=-1),[{0:nmbrivl]) 
title(’Contour Map of Approximation Coefficient Energy Distribution’) 


xlabel([{’Sample number '"n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number (= -Resolution Level) ’) 
pause 


Genvour(ad,10,(0:clmnsize-1), (i2:nmbrivl)) 
title(’Contour Map of Detail Coefficient Energy Distribution’) 


xlabel([{’Sample number "n" (Time = n * ‘,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number (= -Resolution Level) ’) 
pause 

end 


strtsamp = input(’Sample start point? ‘’)+1; 
endsamp = input(’Sample endpoint? ‘’)+1; 
highres = input(’Highest resolution level (least negative)? ’); 
lowres = input(’Lowest resolution level (most negative)? ’); 
endrow = nmbrlivl+l+highres; 
strtrow = nmbrlvl+l+lowres; 
else 
zoom = 0; 
end 
cig 
end 
eic 
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APPENDIX E 


PROGRAM FOR MULTIPHASE ANALYSIS WITH ANY WAVELET 


(hs: 


This program will decompose the input waveform with a Haar wavelet, a 
Daubechies compactly supported wavelet of user-defined order, or with a 
wavelet provided by the user. 

A "sliding" data window is used to maximize signal analysis features 
of the decomposition by negating the effects of random signal time of 
arrival. Approximation and detail coefficients may therefore be 
computed at all resolution levels for every "n" (sample number). 
Because of finite input data length, however, some of the coefficients 
will suffer from "edge effects", i.e., the sample values beyond the 
end of the data sequence are treated as "O"s. 

Output graphs consist of "3~-D" MATLAB mesh plots and "waterfall" 
contour maps of the approximation and detail coefficients, in addition 
to plots of the coefficients at each level (if desired). 


AP AP AY GO HO AP ah? Clo AP GP AP oP ol? 
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Input the "h" coefficients 


BO dO 


elc 
b=[' '}: 
disp(b) 
D1 = [{’ Would you like to use: : 
: 1. your own set of decomposition coefficients, or ” 
4 2. those associated with the Haar wavelet, or : 
. 3. those derived from Daubechies group of wavelets?’]}; 
disp(D2) 
Poeehe— input(’Answer 1, 2, or 3: “*); 
disp(b) 
disp(b) 
if hpick == 1, 
Gisp(’ Note: Sum of coefficients must be normalized to equal 1.7’) 
disp(b) 


hcoeff=input(’What is the name of your decomposition coefficient vector? ‘); 
choice = 1; 

elseif hpick == 2, 
meeerr = (0.5 0.5); 
choice = 1; 

else 
choice = input(’What is the desired wavelet order (2-10 are available)? ’); 
{hcoeff} = daubdata(choice) ; 

end 


hlength = length (hcoeff) ; 


if hlength == 
disp(b) 
disp(’ERROR: Wavelet order selected is outside allowable range.’) 
break 

end 


% Determine if user wants to see plots of the basis functions 
picture=input(/’Do you want to see the Scaling Function/Wavelet (Y or N)? ’,’s'); 
if picture == ‘yY’ 


basisplt 
end 
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& Get the input data vector 


cle 
€Aisp(b) 
Q1 = [’Do you wish to use: J 
’ 1. your own MATLAB formatted row vector, or v 
’ 2. a program generated vector from the following menu?’ 
- Sine wave 4 
: - Pseudo-noise (PN) sequence J 
4 - Sine wave modulated by PN sequence | 
disp(Q1) 


Pick = input(’Answer 1 or 2-3.) 


if Pick ==1, 
% Read in user’s input vector 
Wavedata = input(’What is the name of your input vector? ’); 


sampfreq input(’What was the sampling frequency (Hz)? ’); 
Tsample = 1/sampfreq; 

else 
[Wavedata,Tsample] = wvinput(Pick); 

end 
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hcoeff = sqrt(2)*hcoeff; 


gcoeff fFlipir (neeerr): 
for 97° = 2322hnlengrn 

gcoeff(j}) = -gcoeff(j); 
end 


%* Decomposition algorithm 


datlingth = length(Wavedata) ; 


numrows = ceil(log(datlngth)/log(2)); % determine number of resolution levels 

Lastpts = zeros(1,numrows+1); 

for k = 1:numrows % determine number of points required for each level 
Lastpts(numrows-kt+1) = datlngth+(2*(k)-1)*(hlength-1); 

end 


Lastpts(numrows+1) = datlngth; numpts = Lastpts(1); 
dad = zeros(numrows,numpts); c = zeros(numrows+1,numpts) ; 
c(numrowst1,1i:datlingth) = Wavedata; cworkvct = Wavedata; 
ctr = numrows; 
while ctr >= 1, 

n = numrows-ctr; 

shift = 2*n; oldendpt = Lastpts(ctr+l); newendpt = Lastpts(ctr); 

dnewrow = zeros(1,newendpt); cnewrow = zeros(1,newendpt) ; 

for k = O:hlength-1 coefficient vectors are calculated for each 
resolution level by adding shifted, weighted 
versions of the next higher level "c" vector 

(same effect as direct convolution) 


of oP a oh 


shiftl = k*shift; 
cnewrow(1l+shifti:oldendpt+shift1l) = cnewrow(1+shifti: oldendpt+shifti)+.. 
hcoeff (hlength-k) *cworkvct; 
anewrow(1+shiftl:oldendpt+shift1) = dnewrow(1+shifti:oldendpt+shiftl)+.. 
gcoeff(hlength-k) *cworkvect; 
end 
d(ctr,l:newendpt) = dnewrow; 
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c(ctr,l:newendpt) = cnewrow; 
cworkvct = cnewrow; 
cer e= CLr-1; 

end 


cic 
disp(b) 
D3 = {’ Do you want to see: : 
1. separate plots of the coefficients at each Resolution Level’ 
i in addition to the "3-D" and contour plots, or F 
4 2. Only the “3=D" and contour plots? sa 
disp(D3) 


pickplot = input(’Answer 1 or 2: '); 


% Plot "c" and "a" coefficients by resolution level, if desired. 


if pickplot == 1, 
indxtoO = [O:datlngth-1); 
plotmin = 1.2*min(Wavedata); plotmax = 1.2*max(Wavedata) ; 
if ((plotmin==0) &(plotmax==0)), plotmax = 0.5; end 
if plotmin > 0, plotmin = 0; end 
if plotmax < 0, plotmax = 0; end 
vex {O,datlngth-1,plotmin, plotmax]); 
axis(v); 
plot(indxto0,Wavedata,’*’) 
title(’Approximation Coefficients at resolution 0’) 
xlabel({’Sample number "n" (Time = n * ‘’,num2str(Tsample),’ sec)’])) 
pause 


for k = numrows:-1:1 
n = numrows-k+1; 
endpt = Lastpts(k); 
indxto0O = [0:endpt-1); 
Sroemine = 1.2*min(min(e(k,:)),min(d(k,:)) 
pPlotmax = 1.2*max(max(c(k,:)) ,max(d(k,:)) 
1f ((plotmin==0) &(plotmax==0)), plotmax = 
if plotmin > 0, plotmin = 0; end 
if plotmax < 0, plotmax = 0; end 
v = (0,endpt-1,plotmin, plotmax]; 
axis(v); 
Eubplot( 2 ill) plot(inaxtood,¢c(k,l:endpt) ,"*") 
title(({’Approximation Coefficients at resolution level ’,num2str(-n) } 
Xlabel({’Sample number "n" (Number ef points = ’,num2str(endpt) ,*)’ 
suoplot (212) ,plot (indxto0,d(k,1l:endpt) ,’*") 
title({’Detail Coefficients at resolution level ’,num2str(-n) )) 
Xlabel({’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
pause 


) 
)) 
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% Coefficient energies are used as a measure of response. Display mesh and 
% contour plots of coefficient energy with optimal "zoom" capability to aid 
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Signal analysis 


ao 


C= Cee: 
ds —a.- 2; 
strtsamp = 1; endsamp = numpts; strtrow = 1; endrow = numrows; zoom = 1; 
highres = -1; lowres = -numrows; 
while zoom == 1, 
rangea = [~highres-1:-lowres); ranged = [-highres:-lowres]; 


indxtooO = [strtsamp-1:endsamp-1); 
mesh(c(strtrow:endrowt+1,strtsamp:endsamp) ) 
title(’Energy in Approximation Coefficients’ ) 
if choice -= 1, 
Xlabel({’using Daubechies Wavelet of order ’,num2str(choice) )) 
end 
pause 
contour(c(strtrow:endrow+1,strtsamp:endsamp) ,10,indxto0,rangea) 
title(’Contour Map of Approximation Coefficient Energy Distribution’) 


xlabel({’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number (= - Resolution Level) ’) 
pause 


mesh(d(strtrow: endrow,strtsamp:endsamp) ) 

title(’Energy in Detail Coefficients’ ) 

if choice ~= 1, 

xlabel([’using Daubechies Wavelet of order ’,num2str(choice) }) 
end 

pause 

contour (d(strtrow:endrow,strtsamp:endsamp) ,10, indxto0, ranged) 
title(’Contour Map of Detail Coefficient Energy Distribution’) 


xlabel([’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’}) 
ylabel(’Decomposition Number (= - Resolution Level) ’) 
pause 
cle 
showmore = input(’Would you like to "zoom" in on a section (Y or N)? ’,’s7)5; 
if showmore == /‘Y’, 
disp(b) 
Gdisp(’ You will set the display sample number and resolution level limits.’) 
disp(b) 
aisp({’ Sample range O:’ ,num2str(numpes=1)) 
disp([’ Resolution range -1:’,num2str(-numrows) )) 
revu = input(’Need to see the original contour plots again (Y or N)? ’,’s’); 
if revu == ‘Y’, 


contour (c;10,, [C:numpts—-1),1C:numrows)) 
title(’Contour Map of Approximation Coefficient Energy Distribution’) 


xlabel([’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number = - Resolution Level) ’) 
pause 


contour(d,10, [O:numpts-1), [1:numrows] ) 
title(’Contour Map of Detail Coefficient Energy Distribution’) 


xlabel([’Sample number "n" (Time =n * ’,num2str(Tsample),’ sec)’)) 
ylabel(’Decomposition Number (= - Resolution Level) ’) 
pause 

end 


strtsamp = input(’Sample start point? ’)+1; 
endsamp = input(’Sample endpoint? ‘’)+1; 
highres = input(’Highest resolution level (least negative)? ’); 
lowres = input(’Lowest resolution level (most negative)? ’); 
endrow = numrows+l1thighres; 
strtrow = numrowst+1+lowres; 

else 
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APPENDIX F 


PROGRAM FOR GENERATING ITERATIVE PLOTS OF WAVELETS 
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This program creates and plots the desired iterative approximation 

to the scaling function and the associated wavelet as determined by 
the input "h" vector. The construction is based on the "graphical" 
recursion method. 


Of AO oO oP 


hcoefnew = 2*hcoeff; 
m = length(hcoefnew) ; 


% The "g" vector is determined from the "h" vector. 


gcoefnew = fliplr(hcoefnew) ; 
for jy = 2:2:m 
gcoefnew(j) = -gcoefnew(j); 
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numbrits = input(’How many iterations shall we run for the approximation? '); 


Size = 1; 
hpast = ones(1,m); 
newsize = m; 


for i = l:numbrits 
% Build scaling function approximation 


hnew = zeros(1,newsize); 
for k = 0:size-1 

hnew (2*k+1:2*k+m)=hnew (2*k+1: 2*k+m) +hpast (k+1) *hcoefnew; 
end 


¢ Get wavelet from scaling function approximation and also build time vector 


wv = zeros(1,2*newsize) ;wvlet = zeros(1,newsize);timevctr = zeros(1,newsize); 
shift = newsize/(m-1); 
Bork = O:m-1 
shiftl = round(k*shift+l); 
shift2 = round(newsize+k*shift) ; 
wv(shifti:shift2) = wv(shiftl:shift2)+gcoefnew(k+1) *hnew; 
end 
for k = l:newsize 
wvlet(k) = wv(2*k-1); 
timevctr(k) = k-1; 


end 

interval = (1/2)%*i; 

timevctr = interval*timevctr; % scale time axis 
S = linterval*newsize; 


* Plot basis functions and check calculations with areas and inner products 


if i <= 5, 

plot (timevctr,hnew,'+’,timevctr,wvlet,’*’) 

xlabel({‘Support = ’,num2str(s),’ (Scaling Function: ++++, Wavelet: ****)‘}) 
else 
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plot(timevetr, new, =“, cimevccr, wilet aaa 


xlabel((’Support = ’,num2str(s),’ (Scaling Function: line, Wavelet: dash) ’}) 
end 

titie({‘’Approximations for Iteration number “{numZstr(ae)) 

pause 


hpast = hnew; 
size = newsize; 
newsize = 2*size+m-2; 


phisum = sum(hpast); 


wvysum = sum(wvlet) ; 

sp(i) = phisum*interval; % find area under scaling function 

sw(i) = wvsum*interval; % find area under wavelet 

ip(i)= hpast*wvlet’; % find scaling function/wavelet inner product 


end 


% Output calculation checks 


disp(’ ’) 

aispit’ Area under Area under Inner’) 

disp (’ Scaling Function Wavelet Product’) 

dispt’*) 

fOr 2 = Llenumbrics 

disp¢[‘ f NuUmMZstr (spi) ) / numZ2str(swi));~ ‘ num2str(ip(1)9me 
end 

Gisp(’ ’) 

disp(’<Return>’) 

pause 


clear hcoefnew gcoefnew hpast timevctr wv 
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APPENDIX G 


FUNCTION FOR CONSTRUCTION OF INPUT DATA 
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function [Data,Tsample]) = wvinput(x) 


This function supports the various multiresolution programs by building 
the input test vector desired by the user from the following choices: 
Sine wave, Pseudo-noise sequence, and Sine wave modulated by Pseudo- 
noise sequence with optional additive Gaussian noise. 


of ad? of oP 
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oo 


+ Determine desired input signal type. 
Pigalle 
Q2 = [’ Select the desired type of waveform: f 
J 1. Sine wave f 
: 2. Pseudo-noise (PN) sequence 4 
: 3. Sine wave modulated by PN sequence’); 
ais (5) 
disp (Q2) 


Makev= input (‘Answer 1, 2 Gr 375"); 


% Generate the desired sinusoidal signal. 


if Make ~= 2, 


eic 

% Determine desired sine wave characteristics. 

N2 = [{’ You will determine the following Sine wave parameters: ’ 
- frequency 2 
, - phase ‘ 
a - amplitude ‘ 
- number of samples , 
, - sampling frequency or data length gon ie 

disp(b) 

disp(N2) 

fc = input(’Sinewave frequency ("fc") in Hertz? ’); 


theta = input(’Phase (degrees)? ’); 
A = input(’Amplitude? ’); 
N = input(’Number of Samples ("N", with N = a power of 2)? ’); 
Q3 = [’ Do you wish to set: : 
‘4 1. the sampling frequency, or’ 


: 2. the data length? le 
disp(b) 
disp (Q3) 
SS = anpue(’Answer 1 or 2. 4)> 
1f $5 == 1, 
N3 = [{’ Note: Sampling frequency (fs) must be selected so that 
’ fs/fco >= 2; Data length will be set to allow "N" samples.’]); 
disp(b) 
disp(N3) 
fs2fc = input(’Desired sampling-to-signal frequency ratio? ’); 


GOnSt y= 2*pi7is2tre; 
fs = fs2fc*fc; 


else 
N4 = [{’ Note: Sampling frequency will be chosen to give "N" samples’ 
’ so observation time must be <= "N'"/(2*"fco"). le 
Gisp(b) 
disp(N4) 
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Tobserv = input(’Desired observation time (seconds)? ‘’); 
fs = N/Tobserv; 
const = 2*pi*fc/fs; 
Tsample = 1/fs; 
end 
Tsample = 1/fs; 


$ Construct sine wave 
theta = theta*pi/i8s0; 
HOE kK = 1:N 
Data(k) = sin(const*(k-1)¢theta) ; 
end 
Data = A*Data; 


if Make == 3, % used if PN modulation is desired 
Sinedata = Data; 
Gisp(b) 
Gisp(b) 
N4pt5 = [’ You will now determine the PN sequence characteristics.’ ]; 
disp(N4pt5) 

end 

end 


This section generates the desired pseudo-noise signal by building the 
appropriate feedback shift register (FSR). 


oO” of? 


if Make -= 1, 
eric 
% Determine desired PN sequence characteristics. 
N5 = {’ Note: PN sequence = {+1 or -1,...,2*m - 1) 


Feedeack, shift register initial state is {1,1,...,1)']; 
disp (b) 


disp(N5) 
N6 = {’ You will determine the following parameters: ’ 
- number of stages ; 
: - number and position of taps : 
U - chip rate 
y - time delay : 
J - number of samples 
: - sampling frequency or data length le 
disp(b) 
aisp(N6) 


mstages = input(’Desired number (m) of stages (2-10)? ’); 
numbrtap = input(’Number of taps? ’); 
for k = l:numbrtap 
Tap(k)=input([{’What is the position of Tap number ’,num2str(k),’? ‘’)); 
end 
chiprate = input(’Chip rate (chips/sec)? ’); 
delay = input(’Sequence delay in chips (positive real number)? ’); 
if Make == 2, 


If the PN signal alone is desired the user must input the number 
of samples and sampling rate; otherwise, they are dictated by 
the choices made for the sinusoid. 


o° co oO 


me— Sort(2); 
N = input(’Number of samples? ’); 
Q4 = [’ Do you wish to set: ' 
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i 1. the sampling frequency, or’ 


7 2. the data length? ae 

disp (b) 
disp (Q4) 
P6 = input(’Answer 1 or 2. ’); 
1f P6é == 1, 

N7 = [’ Note: Sampling frequency ("fs") must be selected so that ’ 

’ fs/chiprate >= 2; data length will be set for "N" samples.’); 
disp(b) 
disp(N7) 


fs2cr = input(’Desired sampling frequency to chip rate ratio? ’); 
Tsample = 1/(chiprate*fs2cr); 


else 


N8 = [’ Note: Sampling frequency will be chosen to give "N" samples’ 
’ so observation time must be <= "N'"/(2*Chip rate). 


disp(b) 
disp (N8) 


Tobserv = input(’Desired observation time (seconds)? ’); 


fs = N/Tobserv; 
fs2cr = fs/chiprate; 
Tsample = 1/fs; 
end 
else 
disp(b) 


N9 = [’ Same number of samples (’,num2str(N),’) and sampling ’); 
N10 = [’ frequency (’,num2str(fs),’ Hz) that were selected for’); 
N11 = [’ the sine wave will be used. <Return>’]; 


disp(N9) 

aisp(N10) 

disp(N11) 

pause 

fs2cr = fs/chiprate; 
end 


% Generate base PN sequence 
FSR = ones(1,mstages) ; 

Lb = 2°mMstages = 1; 

cklength = N/(fs2cr*L); 
repeat = ceil(cklength) ; 

V = zeros(l1, (repeattl) *L) ; 
forsKk = 12L 


V(k) = FSR(mstages) ; 
sigma ; 
for m 


‘ 


l:numbrtap 


Sigma = sigma+FSR(Tap(m) 


end 


% Determine how many periods are necessary 


% build one period of the sequence 


 modulo-2 tap adder 
Ve 


FSR(2:mstages) = FSR(1l:mstages-1); 


FSR(1) = rem(sigma,2); 
end 


© 


for n = l:repeat 


% Ensure enough periods are available for the desired number of samples 


NV (n* bi tnt) *b) ay Cee, 


end 
delay = rem(delay,L); 
for n- = 12:N 

Data(n 


) = V(f£ix((n-1)/fs2cr + delay)+41) ; 
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as 


if Data(n) == 0, Data(n) = -1; end % make sequence bipolar 
end 


% If selected, modulate the sine wave with the PN sequence. 


if Make == 
Data = Sinedata.*Data; 
end 


$ Plot the constructed signal and add noise, if desired. 


indxtoo = [0:N-1]}; 

v = (0,N-1,-1.2*max (Data) ,1.2*max (Data) }; 

axis(v); 

plot (indxto0,Data) ,title(’MATLAB Plot of Input Signal Vector’) 
xlabel(’Sample number’) ,ylabel(’Volts’) 

pause 

axis; 


Src 

disp(b) 

disp(’ White Gaussian Noise is generated by the MATLAB "random" function. ’) 
noise = input(’Do you want noise added to your signal (Y or N)? ’,’s’); 

if noise == ‘Y’, 


% SNR based on "continous" signal energy and the Gaussian noise variance 


SNRGB = input(’/What is the desired Signal-to-Noise Ratio (dB)? ’); 
SNR = 10*(SNRdGB/10); 
Noisepwr = ((A%2)/2)/SNR; 
rand(’seed’,0); rand(’normal’); 
noisvctr = Noisepwr*rand(1,N); 
Data = Datat+noisvctr; 
we O, N-1,-1.2*max (Data) ,1.2*max (Data) j; 
axis(v); 
plot (indxto0, Data) ,title(’MATLAB Plot of Input Data Vector’) 
Xlabel(’Sample number’) ,ylabel(’Volts’ ) 
pause 
axis; 
end 


return 


eat 


APPENDIX H 


FUNCTION FOR DAUBECHIES' H-COEFFICIENTS 
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function{hcoeff] = 


daubdata (x) 


% This function returns the proper "h" coefficients for the specified order 
% Daubechies wavelet. The input value of "x" is the specified order. 


hcoeff 
elseif x 
hcoeff 


elseif x 
hcoeff 


elseif x 
hcoeff 


elseif x 
hcoeff 


elseif x 
hcoeff 


elseif x 
hcoeff 


elseif x 
hcoeff 


elseif x 
hcoeff 


else 
hcoeff 
break 
end 


hcoeff = 


Gerurn 


emae2962914145,0-836516303738;0.224143868042;-0.129409522551)] ; 
| 


(03332670552950;0.806891509311;0. 459877502118; -0.135011020010; 

SWePUCo4a 27 SeGe nOel so ezocg1ceZ |); 

== 4 

eepo.2503776135309;,0.71484657055370.630880767930;-0.027983769417; 
[Oo Uso too, 0.05004 17381536;0.032883011667 ;-0.020597401785)}; 

—— > 

SeeOe. bout 2569797450. 603e2oo797;0.72430852E8438°0.138428145901;; 
—-0.242294887066;-0.032244869585;:0.077571493840;-0.006241490213; 
UO oo Ove, 0. 002355725265); 

== 6 

=eOe da 5407455 50;02-494623890398 -0.751133908021;0.315250351709; 
POmeezocuaorsdGo, = O.)2975656675607 ;0.097501605587 ;0O.027522865530; 
museowooezuss slo 7O 000553642201 70.,004777257511;-0.001077301085] ; 

—— ne), 

= {(0.077852054085;0.396539319482;0.729132090846;0.469782287405; 
-0.143906003929;-0.224036184994 :;0.071309219267;0.080612609151; 
wo soOgZe2oOsI5,—-0. 0165749041631; 0.,012550998556;0.000429577973; 
a0 001601640704 ;O20003537123800] ? 

=—— 

si eUoaeaioe4224a5 70231228 71590914;0.675630736297;0.585354683654; 
=0. 015829105256; =-0. 28401554 2962;0.000472484574;0.128747426620; 
=O. 017369301002 ;—-0.044088253931;0-013981027917;0.008746094047; 
=-0.004870352993;-0.000391740373;0.000675449406;-0.000117476784); 

== 9 

Siow eeo07 7247326670. 243634674615;0.6042623123690;0.657288078051; 
Ores oto 7 262020 77022952717 5783279;—-0.096840783223;70.148540749338; 
Gewisw7 25051479, —-0, 067652629061 ;0.000250947115;0.022361662124; 
woe 004272220475e;, —-0-004281503682;7;0.001847646883 ;0.000230385764; 
a0. 0002517563789; 0. 000039347320); 

= {[0.026670057901;0.188176800078;0.527201188932;0.688459039454; 

Opec 7254 061, -0.249646424327;-0.195946274377;0.127369340336; 

0.093057364604 ;-0.071394147166;-0.029457536822;0.033212674059; 

P7002 SCs n 5 5o7;,—-0.010733175483 ;0.001395351747;0.001992405295; 

=O7000GGSG00605,-0.000116466855;0.000093588670;-0.000013264203); 


sible 


heoetf’ /sart(2) ; % normalize the "h" vector 
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